function U = quasistatic_MT_greensfunc(mdl,obs,mu,nu)
%coded by josh crozier 6/12/2018
%quasi-static moment tensor greens function
%outputs displacement U
%same input format as disloc3d okada code, except:
    %mdl = [sourcedepth sourceeast sourcenorth EE NN ZZ EN EZ NZ]
%singleforce functions from Segall EQ and Volcano deformation pg 77-78

%force couples (MT components)
M11 = mdl(4);
M22 = mdl(5);
M33 = mdl(6);
M12 = mdl(7);
M13 = mdl(8);
M23 = mdl(9);

%source loc
e1 = mdl(2);
e2 = mdl(3);
e3 = mdl(1);

d = 10^(-2); %force couple offset, just needs to be small (<1m fine)

%single forces offset by +-0.5 in east dir
mdltemp = [e3, e1+d/2, e2, 1, 0, 0];
gpos1_1 = singleforce_greensfunc(mdltemp,obs,mu,nu);
mdltemp = [e3, e1-d/2, e2, 1, 0, 0];
gneg1_1 = singleforce_greensfunc(mdltemp,obs,mu,nu);

mdltemp = [e3, e1+d/2, e2, 0, 1, 0];
gpos1_2 = singleforce_greensfunc(mdltemp,obs,mu,nu);
mdltemp = [e3, e1-d/2, e2, 0, 1, 0];
gneg1_2 = singleforce_greensfunc(mdltemp,obs,mu,nu);

mdltemp = [e3, e1+d/2, e2, 0, 0, 1];
gpos1_3 = singleforce_greensfunc(mdltemp,obs,mu,nu);
mdltemp = [e3, e1-d/2, e2, 0, 0, 1];
gneg1_3 = singleforce_greensfunc(mdltemp,obs,mu,nu);

%single forces offset by +-0.5 in north dir
mdltemp = [e3, e1, e2+d/2, 1, 0, 0];
gpos2_1 = singleforce_greensfunc(mdltemp,obs,mu,nu);
mdltemp = [e3, e1, e2-d/2, 1, 0, 0];
gneg2_1 = singleforce_greensfunc(mdltemp,obs,mu,nu);

mdltemp = [e3, e1, e2+d/2, 0, 1, 0];
gpos2_2 = singleforce_greensfunc(mdltemp,obs,mu,nu);
mdltemp = [e3, e1, e2-d/2, 0, 1, 0];
gneg2_2 = singleforce_greensfunc(mdltemp,obs,mu,nu);

mdltemp = [e3, e1, e2+d/2, 0, 0, 1];
gpos2_3 = singleforce_greensfunc(mdltemp,obs,mu,nu);
mdltemp = [e3, e1, e2-d/2, 0, 0, 1];
gneg2_3 = singleforce_greensfunc(mdltemp,obs,mu,nu);

%single forces offset by +-0.5 in z dir
mdltemp = [e3-d/2, e1, e2, 1, 0, 0];
gpos3_1 = singleforce_greensfunc(mdltemp,obs,mu,nu);
mdltemp = [e3+d/2, e1, e2, 1, 0, 0];
gneg3_1 = singleforce_greensfunc(mdltemp,obs,mu,nu);

mdltemp = [e3-d/2, e1, e2, 0, 1, 0];
gpos3_2 = singleforce_greensfunc(mdltemp,obs,mu,nu);
mdltemp = [e3+d/2, e1, e2, 0, 1, 0];
gneg3_2 = singleforce_greensfunc(mdltemp,obs,mu,nu);

mdltemp = [e3-d/2, e1, e2, 0, 0, 1];
gpos3_3 = singleforce_greensfunc(mdltemp,obs,mu,nu);
mdltemp = [e3+d/2, e1, e2, 0, 0, 1];
gneg3_3 = singleforce_greensfunc(mdltemp,obs,mu,nu);

%MT greens functions
g11 = (-gneg1_1 + gpos1_1)/d;
g22 = (-gneg2_2 + gpos2_2)/d;
g33 = (-gneg3_3 + gpos3_3)/d;
g12 = (-gneg2_1 + gpos2_1 -gneg1_2 + gpos1_2)/d;
g13 = (-gneg3_1 + gpos3_1 -gneg1_3 + gpos1_3)/d;
g23 = (-gneg3_2 + gpos3_2 -gneg2_3 + gpos2_3)/d;

%displacements
U= M11*g11 + M22*g22 + M33*g33 + M12*g12 + M13*g13 + M23*g23;
end

